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ABSTRACT 

We model the dynamical evolution of primordial black holes (BHs) in dense star clusters using a simplified 
treatment of stellar dynamics in which the BHs are assumed to remain concentrated in an inner core, completely 
decoupled from the background stars. Dynamical interactions involving BH binaries are computed exactly and 
are generated according to a Monte Carlo prescription. Recoil and ejections lead to complete evaporation of the 
BH core on a timescale ~ 10 9 yr for typical globular cluster parameters. Orbital decay driven by gravitational 
radiation can make binaries merge and, in some cases, successive mergers can lead to significant BH growth. 
Our highly simplified treatment of the cluster dynamics allows us to study a large number of models and to 
compute statistical distributions of outcomes, such as the probability of massive BH growth and retention in 
a cluster. We find that, in most models, there is a significant probability (~ 20-80%) of BH growth with 
final masses > 100M Q . In one case, a BH formed with mass 620M q . However, if the typical merger 
recoil speed (due to asymmetric emission of gravitational radiation) significantly exceeds the cluster escape 
speed, no growth ever occurs. Independent of the recoil speed, we find that BH-BH mergers enhanced by 
dynamical interactions in cluster cores present an important source of gravitational waves for ground-based 
laser interferometers. Under optimistic conditions, the total rate of detections by Advanced LIGO, could be as 
high as a few tens of events per year from inspiraling BHs from clusters. 

Subject headings: galaxies: star clusters — globular clusters: kinematics and dynamics — black hole physics — 
gravitational waves 



1. INTRODUCTION 

1.1. Astrophysical Motivation 

Many observations of globular clusters suggest the possi- 
ble existence of intermediate-mass black holes (IMBHs) with 
masses ~ 10 3 - 10 4 M Q in the centers of some cluster cores. 
The predicted masses of these IMBHs agree well with a sim- 
ple extrapolation of the M - a (mass - velocity dispersion) 
relation for galactic nuclei JGebhardtetalJ l2000). Observa- 
tions and dynamical modeling of the globular clusters Ml 5 
in the Milky Way and Gl in M31 appe ar to be con sistent 
with such a central massive BH dGerssen et al.l l2002, 2003; 
IGebhardt., Rich., ffi; Ho||20 Q2 | ). However, TV-body simulations 
by Baumgardt et al. (2003a b) suggest that the observations 
of M15 and Gl, and, in general, the properties of all core- 
collapsed clusters, could be explained equally well by the 
presence of man y compact stars near the center without a 
massi ve BH (cf. Ivan der Mareill2004t IGebhardt. Rich. & Hoi 
120051) . On the other hand, these A^-body simulations also 
suggest that many, perhaps most, non-core-collapsed clus- 
ters (representing about 80% of globular clusters i n the Milky 
Way) could contain a central IMBH (Baumgardt et al. 2004b, 
2005). 

Although the origin of IMBHs is not directly constrained 
by any observations, one possibility that has received consid- 
erable attention is the growth of a very massive object through 
successive collisio ns and mergers of ordinary massive stars 
in the cluster core (Colgate 1967; Ebisuzaki et al. 200 1 |. Re- 
cent numer ical studies bv iPor tegies Zwart & McMillan 
(20 021), IGiirkan. Freitag & Rasiol (120041) . and 
Freita g~Gurkan. & Rasiol il2005l) demonstrate that mass 
segregation and core collapse could proceed so quickly that 
there is rapid, runaway growth through stellar collisions 
before the most massive stars have evolved to supernovae 
(in < 3Myr). Other numerical evidence comes from the 



iV-body simulations of Porte gies Zwart et"afl (|2£)04), which 
demonstrate such growth in some young clusters. 

An important alternative, which we study in this paper, is 
the growth of an IMBH thro ugh successive mergers of stel- 
lar ma ss (~ 10M Q ) BHs (Ouinlan & Shapiro 19891 lLeell995l 
1200 ll) . In this case, the massive stars formed initially in the 
cluster (with masses > 20M©) must have avoided physical 
collisions and mergers and instead were able to complete their 
normal stellar evolution and produce BHs. Unlike massive 
stars, BHs have negligible cross sections for direct collision, 
so BH mergers can only occur through gravitational wave 
(GW) emission in binaries, possibl y enhanced by dy namical 
interactio ns in the cluster (see Miller & Colbert 2004 for a re- 
view and Sll.2> . 

Even if it does not lead to growth and IMBH formation, 
the dynamical evolution of BHs in dense star clusters could 
also produce large numbers of tight BH-BH binaries that will 
merge through GW emission (possibly outside the cluster fol- 
lowing dynamical ejection). These merging BH-BH binaries 
are likely candidates for direct GW detection and could even 
dominate the detection rates for ground-based laser interfe r- 
ometers such as LIGO ( Porte gies Zwart & McMillanl2 000). 

1.2. BH Formation and Segregation 

In a globular cluster, it can be expected that a frac- 
tion — 10" 6 to 10" 4 of the TV ~ 10 6 initial stars will 
become stellar-mass BHs v ia normal stellar evolution 
(Sigurdsson & Hernquist 1993). Assuming that all stars with 
initial mass greater than 20M Q become BHs, the most re- 
cent K roupa initial mass function (IMF) ( Kroup a & WeidneJ 
2003) gives a slightly higher initial BH fraction of jVbh ~ 
1 .5 x 10" 3 jV, where we use m m ; n = .08M Q and m max = 150M Q 
as the minimum and maximum stellar masses. When we 
scale this to the total mass of the cluster, M c \, we find A^h ~ 



2 



O'Leary et al. 



3x10 3 (M c i/M Q ). All of these BHs should have formed be- 
fore ~ 10 Myr, w ith the most massive BHs forming at around 
3 Myr ( Schalle retaiT!99l . 

The radial distribution of these BHs in the cluster when 
they form is not known, but it is reasonable to assume that 
they should be much more centrally concentrated than the re- 
maining main-sequence stars (MSs). This is for three rea- 
sons: (1) we expect significa nt mass se gregation of the ini- 
tial, higher-mass progenitors (Freitag et al. 2005); (2) there 
may be preferential formation of massive stars near the clus - 
ter center (see, e.g., Mur ray & Linll996tlBonnell et all 2001 ): 
(3) BH birth kicks are not expected to be large enough to dis- 
place the BHs into the cluster halo (or eject them from th e 
cluster; see White & van Paradiisll996llWillems et all2005l) . 
Even if the BHs were initially distributed throughout the clus- 
ter, mass segregation would still be able to concentrate them 
into a central sub-cluster on a relatively short timescale. In- 
deed, a BH of mass Mrh near the half-mass radius will be 
brought into the cluster core on the timescale 



Mr 



-trh, 



(1) 



where t± is the relaxation time at the half-mass radius, and 
(m) is the average stellar mass (Fregeau et al. 2002). Consid- 
ering a typical cluster with t± ~ 1 Gyr, we see that a sub- 
cluster of BHs should still assemble near the center after at 
most — 100 Myr. 

After a time very short compared to the overall dynamical 
evolution timescale of the globular cluster, the BHs should 
then form a self-gravitating subsystem within the core, which 
is dynamically decoupled from the rest of the cluster. This 
decoupling is sometimes referred to as Spitzer's "mass strati- 
fication instability" (Spitzer 1969). The physics of this insta- 
bility is by now very well understood, both for simple two- 
component ^y_stemsand for clusters with a b road mass func- 
tion ilWatterset all2000t IGurkan et all2"0 04 ) . 

The dynamical evolution of the BH subsystem proceeds 
on a much shorter timescale since its relaxation time is now 
~ Nbh/N times smaller than for the parent cluster. Interac- 
tions involving hardening of BH binaries will eventually lead 
to the ejection of BHs from the cluster, until there are so few 
left that the BH subsystem recouples dynamically (returns 
to "thermal equilibrium") with the other cluster stars, and 
the interaction rate between BHs and cluster stars becomes 
comparable to the BH- BH interact ion rate. For simple two- 
component clusters, Spitzer ( 1969) derived through analytic 
methods the condition necessary to reach energy equiparti- 
tion: 



where Mi < Mi are the total masses of the two compo- 
nents, and m? > m\ the ma ss of each individual object. 
Watte r~Toshi. & Rasiol (|2000) used an empirical approach to 
find the more accurate condition 



(^r<0.32. 



(3) 



For a cluster with total mass Mi = 10 6 Mq, BH mass W2 = 
15 M©, and average star mass mi = 1 M Q , the cluster can be in 
equipartition, according to equation Q, when there are < 30 
BHs in the cluster, i.e., significantly fewer than the number 
expected to form from the IMF. The minimum number of BHs 



required to decouple from the cluster would li kely be even le ss 
if their mass spectrum was considered (Giir kan et all2 004). 

For most of its subsequent dynamical evolution, we expect 
the BH sub-cluster to remain largely free of other kinds of 
stars, even MSs with comparable masses. Indeed, for BHs 
to have formed, the most massive stars in the cluster must 
have avoided runaway collisions. T his requires an initial half- 
mass relaxation time t± > 30Myr ( IGurkan et alJl2004l) . On 
the other hand, by the time the MS turnoff has decreased 
to IOMq, driving core collapse to concentrate the remain- 
ing IOMq MSs into the core (where they could then interact 
with the B Hs) befor e they evolve would require t± < 20 Myr 
iGurkan et alJ 120041 see especially their Fig. 10). There- 
fore we see that the parameter space for clusters where both 
BHs and massive MSs would segregate and decouple simul- 
taneously is probably very small, or nonexistent. This will 
allow us to concentrate on "pure" BH systems in our nu- 
merical simulations (§2). This picture is not altered signifi- 
cantly by the presence of binaries. During mass segregation, 
as the BHs start concentrating into the denser cluster core, 
they will likely interact with each other to form BH-BH bi- 
naries, if they were not already in binaries. BH-MS binaries 
will quickly undergo three-body or four-body exchange inter- 
acti ons that replace t he lighter MS companion by a heavier 
BH ( Sigur dsson & Phinnevl 19931) . Therefore, as the BH sub- 
cluster begins to dynamically decouple from the rest of the 
cluster, a considerable number of hard BH-BH binaries are 
expected to remain. This has been demonstrated with direct 
iV-body simulations by Portegies Zwart & McMillan (2000), 
where they found that almost all BH-binaries ejected were 
BH-BH. 

1.3. Previous Studies 

Motivated by the absence of BH X- ray bi- 
naries, Kulkarni, Hut. & McMiiiar] Jl993h and 
Sigurdsson & Herna uisa d!993l) discussed the evolution 
and fate of 10M Q BHs in globular clusters using simple an- 
alytic considerations. Their conclusions can be summarized 
as follows. Through dynamical interactions, hard BH-BH bi- 
naries tend to be hardened further (whereas soft binaries tend 
to be disrupted; lHeggid ll975). Eventually, as the BH-BH 
binaries harden, the recoil produced by interactions becomes 
so great that the binaries can be ejected from the cluster. The 
timescale for merger by GW emission usually remains longer 
than the interaction time in the BH cluster so that hardened 
binaries are almost always ejected. Eventually the number of 
BHs is depleted, and no more than a few BHs would remain 
in the cluster core. This would imply that BH growth through 
successive mergers in the cluster cannot occur. Even if ~ 1 
BH remained at the center of every globular cluster today, 
it is unlikely that this would ever become d etectable as an 
X-ray binary (Kalogera, King, & Rasio 2004). 

In order to bett er under stand the complex int eractio ns of 
BHs in clusters, iPortegies Zwart & McMillanl feOOOl) per- 
formed direct A^-body simulations of systems with N = 2048 
and N = 4096 total stars, including a small fraction (~ 1%) 
of equal-mass "BHs" 10 times more massive than the other 
stars. They found that ~ 90% of the BHs were ejected from 
the cluster after 4-10 relaxation times of the cluster (less 
than a few Gyr for most clusters), including ~ 30% in BH- 
BH binaries. Simila r Af-body s imulations were most recently 
performed by Merritt et alJ (12004). who studied the formation 
of a core of BHs in a two-component cluster, with N = 10 4 . 
iMerritt et alJ ( | 2004) found that the BHs completely segregate 
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to the core within ~ lOOMyr (ignoring initial mass segrega- 
tion of the higher mass progenitors). These results are consis- 
tent with our qualitative understanding of the interactions of 
BHs, and strongly supports the assumptions we will make in 

§12 

Several studies proposed scenarios that could help BH- 
BH binaries merge inside clusters, opening the possibility of 
BH growth and IMBH formation through successive mergers. 
Miller & Hamilton ( 2002b) suggested that one larger seed BH 
may help overcome the Newtonian recoil, since the mass of 
the binary would be too large to have a recoil velocity greater 
than the escape velocity from the cluste r, thus anchoring it 
in the c ore. It has also b een proposed by M iller & Ham ilton 
(2002aj) and lWenl (12003) that binary-binary interactions may 
have a large influence on BH mergers in the cluster core. 
These authors suggest that binary-binary interactions can pro- 
duce significant numbers of long-lived hierarchical triple sys- 
tems in which the outer BH increases the in ner binary's eccen - 
tricity via Kozai-type secular perturbations llFord et alJ' 2000). 
thereby increasing the merger rate. Because these hierarchi- 
cal triples may be driven to merge before their next interaction 
they should have a higher probability of staying in the cluster, 
a nd this can be a mechanism for retaining merging BHs. 

iGultekin. Miller. & Hamiltonl (12004 (hereafter GMH04) 
were the first to look at the possibility of successive merg- 
ers of BHs in a cluster core. In their simulations, GMH04 
computed successive interactions between a BH-BH binary 
of varying mass ratio with single 10M Q BHs sampled from 
an isotropic background. Between interactions the binary 
was evolved according to general relativity. They concluded 
that GW emission allowed for more mergers than previously 
thought possible, but the number of BHs required to form an 
IMBH would be much greater than believed to exist in a typi- 
cal globular cluster. 

In our simulations, we treat not only binary-single interac- 
tions, but also four-body (binary-binary) interactions. Most 
importantly, we compute dynamical interactions for a realis- 
tic mixture of single and binary BHs self-consistently within 
a "pure BH" cluster core. We implement a treatment for the 
secular evo l ution of triples, including the Kozai mechanism 
dWenll2003UMiller & Hanultonll2002al) . We also include, for 
the first time in any dynamical treat ment, a more realistic IMF 
for the BHs (to be detailed in 12 .2\ . 

Our paper is organized as follows. In 50 we present our 
simulation methods and assumptions, as well as all our initial 
conditions. The main results of our simulations are shown in 
<2] where we look at evolution of all BH sub-clusters, and in 
|4] where we look at merging binary BHs as a source of GWs. 
Finally we conclude our paper in 5j5]with a summary and dis- 
cussion of the implications of our results, and suggestions for 
further studies. 

2. METHODS AND ASSUMPTIONS 

2.1. Numerical Methods 

Ideally, to simulate the evolution of ~ 10 3 BHs in a massive 
cluster, one would like to do a full iV-body simulation of the 
entire cluster, including the BHs and other stars. Or, capital- 
izing on the fact that the BH subsystem dynamically decou- 
ples from the rest of the cluster early on, one could perform 
an TY-body simulation of just the BHs, subject to the cluster 
potential due to the other stars and the effects of dynamical 
friction, which tend to bring BHs that have been kicked out 
of the core — but not out of the cluster — back into the core on 
a short timescale. Although the second scenario is evidently 



computationally much cheaper than the first to treat via di- 
rect TY-body techniques, we decided to adopt an even faster 
technique which, although approximate, includes all the vi- 
tal physics of the problem. This allows us to sample a wider 
range of initial conditions, and make a more thorough map of 
the parameter space of the problem. 

We treat the evolution of a BH subsystem in a background 
cluster subject to binary interactions, using a Monte Carlo 
technique to sample interaction rates and treat mass segrega- 
tion, in conjunction with the small iV-body tool kit Fewbody 
to nu merically integrate binary interactions dlvanova et alJ 
120051) . We assume that the BH sub-cluster has a constant 
density and velocity dispersion throughout its evolution. The 
justification for this assumption is that with any reasonable 
initial binary fraction (> a few percent), the sub-cluster will 
enjoy a long-lived binary-burning phase in whic h its core pa - 
rameters are roughly constant with time (Fregeau et al. 2003). 
The code we use is a modified version of the one presented in 
llvanova et alJ d2005). specially adapted to treat a BH system. 

Each dynamical binary interaction is followed using 
Fewbody, a numerical toolkit for simulating small-TY grav- 
itational dynamics that provides automatic calculation termi- 
nation and classificatio n of o utcomes (for a detailed descrip- 
tion see Fregea uet alJ 120 04). Fewbody numerically inte- 
grates the orbits of small-TY systems, and automatically clas- 
sifies and terminates calculations when an unambiguous out- 
come is reached. Thus it is well-suited for carrying out large 
numbers of binary interactions, which can be quite complex 
and long-lived, and thus must be treated as computationally 
efficiently as possible. 

Star clusters are characterized by a dense central core, sur- 
rounded by a much larger low-density halo. Consistent with 
this structure of a star cluster, we assume that all strong inter- 
actions occur within the core, which has a velocity dispersion 
Ccore (cf. eq. 11 II ). If a product of an interaction has a velocity 
greater than the escape velocity from the core of the cluster, 
v eS c, then it is assumed ejected from the cluster and is removed 
from the simulation. If it is less than v esc but greater than the 
escape velocity from the core into the halo Vhaio, it is placed 
in the halo of the cluster, from where it can later reenter the 
core through dynamical friction with the background stars. 
Dynamical friction is implemented in our code by sampling 
from a Poisson distribution with an average timescale given 
in equation Q with the average s tellar mass, (m) = 1 .0M© 
(see §3.3 from llvanova et all20 05). 

Between interactions, all BH-BH binaries are evol ved ac- 
cordin g to the standard post-Newtonian equations (iPetersI 

da 64 G 3 m\ m\(mi + mj) ( 73 2 37 A 
~dt = ~~5~^ fl 3 (1 _ e 2 )7 /2 [ 1 + 24 e + 96 S J () 

de 304 G 3 mi mi(m\ +ni2)e / 121 2 \ 
dt = ~l5^ a 4 (l-e 2 ) 5 /2 \ + 304 e )-' (5) 

where m\ and m-i are the masses of the two BHs, a is the 
binary semimajor axis, and e is the orbital eccentricity. 

In some simulations, we account for linear momen- 
tum kicks imparted to t he binary due to the asymme- 
try of the GW emission (Fitchet q 119831). Because of the 
large theoretical uncertain ty (Favata, Hugh es. & Holzll2004t 
Blanchet, Ousailah, & Will 2005) in the recoil velocity of 
"major mergers" (i.e., when the mass ratio q = niilm\ > 0.4; 
here we assume m\ > mi), and the smaller but significant un- 
certainty from the spins of the BHs, we opted to neglect spin 
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in determining recoil velocities in our simulations. We deter- 
mine the overall recoil velocity of the merg er remnant, V r ^ 
by using the form of the equation derived by Fitchett ( 1983), 



J max \ Isco 



(6) 



where f(q) = q 2 (\-q)/(\+q)\ f max ss .0179, V is the maxi- 
mum magnitude of re coil, and n sm is the radius of innermost 
stable circular orbit. iFitchettl d!983h found for circular or- 
bits Vo w 1480 km s" 1 , much great er than the escape v elocity 
from any globular cluster, whereas Favata et al. (2004) found 
Vo ~ 10 - 100 km s" 1 . In our simulations we set 



/(«) 

J max 



(7) 



for ease of comparing Vo with other works. We use Vo near or 
slightly above the escape velocity of the cluster, up to 80 km 
s" 1 . The form of equ ation (0 is consistent with the analysis 
of ( Blanchet et al. 2005 ), who found the recoil velocity for the 
merger of two non-spinning BHs to be Vo w 250 ±50 km s" 1 at 
the second post-Newtonian order. We do not account for GW 
recoil in the merger of the inner binaries that are part of hier- 
archical triples, because in the simulations where we include 
G W recoil, m ergers in hierarchical triples are insignificant. 

iLed (119931) shows that for the velocity dispersions ( < 100 
km s" 1 ) and numbers of BHs ( < 10 3 ) expected in the star clus- 
ters we are investigating, the rate of two-body binary forma- 
tion from gravitational bremsstrahlung is much less than that 
of regular (Newtonian) three-body binary formation, whereby 
a binary is formed with the help of a third BH, which takes 
away the excess energy needed to form the bound pair. There- 
fore, for our simulations, we only account for three-body and 
binary-binary (four-body) Newtonian interactions. 

In a dense sub-cluster of BHs, three-body binary formation 
can lead to the formation of a significant number of BH bi- 
naries, and eventu ally can help lead to the disruption of the 
entire sub-cluster. Ilvanova et al.l (12005) calculated the three- 
body binary formation rate for a binary of minimum hardness 



'/n 



(8) 



where fo max is the maximum size of the region in which the 
three objects interact and cr core is the three-dimensional veloc- 
ity dispersion of the core. The final rate per star they found 
for the formation of a binary with hardness r\ > rj m i n is 



nlGHm) 5 



T(rj > ?7min) = 7T- 

<7™ 



- /(mi, m 2 ,m 3 , 77) 



(9) 



where 



a k «2«3 m\ m\ _ 5 

f(mi,m 2 ,m3,ri)= — T - — -±r— -^77 (I+277) 



n c (m) 5 (m) 5 



1 + 



V3 



-V 



-1/2 



mi '«2 



(ml +m2) (m) / a, 



V12 



(10) 



n c is the core density of the BHs, n 2 is the core density of 
objects of mass mi, m is the core density of objects of mass 
m3, V12 is the relative velocity of the first object with respect 
to the second, and V3 is the relative velocity of the third ob- 
ject with respect to the center of mass of the first two objects. 
We follow the exac t treatment of dynamical interactions of 
Ilvanova et al.l (12005) but include three-body binary formation 



with T/min = 1 in a consistent manner (see, in particular, their 
§3.4). We note that a more accurate criterion for the minimum 
binding energy we use should be based on the orbital velocity 
of the lightest member of the formed binary since we are not 
looking at equal-mass clusters (Hills 1990). However, since 
we typically do not have mass ratios above ~ 10, we find our 
criterion to be sufficient. 

Because the code is not yet capable of following the evo- 
lution of triples between interactions, we must break them 
up before the next interaction time-step. In order to de- 
termine how to destroy the triple, we check if the binary 
is likely to merge before its next interaction. As a first 
step, we integrate numerically equations and (|5} (iPetersI 
1964ft. We also scale the timescale of merger according to 
Miller & Hamilton (2002a) by calculating the maximum ec- 
centricity from a first order Kozai mechanism approximation 
without pos t-Newtonian precession by solving equation (8) of 
IWenl (120031) . We then use the smaller merger time of the two 
methods. It is necessary to consider both methods because 
the scaling from Miller & Hamilton ( 2002a) overestimates the 
merging time in the instances when the Kozai mechanism is 
insignificant. In this case, the inner binary is merely perturbed 
and the eccentricity does not fluctuate, therefore the timescale 
of the merger should be the same as for an unperturbed binary. 
If the inner binary is likely to merge before the triple would 
interact with a field BH or BH-BH binary it is immediately 
merged, otherwise the triple is broken-up keeping the inner 
binary. Here, we assume that the outer BH is ejected from the 
triple with a velocity equal to its relative velocity with the cen- 
ter of mass of the triple. We keep the inner binary, but shrink 
the binary separation to conserve the energy of the system. 

2.2. Initial Conditions and Parameters 

We use the results of iBelczvnski. Sadowski. & Rasid 
(2004) (hereafter BSR04), adopting the mass and binary dis- 
tributions of their standard model at ll.OMyr for our initial 
conditions (see, e.g., their Fig. 2, Fig. 4, & Fig. 5). BSR04 
used a population synthesis approach to follow the evolution 
of a large number of massive stars and binaries, as would 
likely form in a massive star cluster. The model we base our 
calculations on has an initial binary fraction /,;„ = 50%, and 
follows the traditional Salpeter IMF for all initial stars with 
mass > 4Mq. The binary fraction we use for our simula- 
tions is, of course, the final binary fraction found in BSR04, 
/bin ~ 14%. Although most BH binaries have MS compan- 
ions, we assume that these binaries will eventually become 
BH-BH binaries through exchange interactions. We create 
BH-BH binaries in their place, and select the companion BH, 
such that the distribution of the mass ratio, q, is uniform 
throughout the range q m \ n < q < 1 (q m \ n is set by the mini- 
mum BH mass in our dis tribution) , consistent wit h observa- 
tions for q > 0.2 (iWoitas. Leinert. & KohTerjl2001l) . We then 
increase the separation of the BH-BH binary assuming the 
binary preserves its binding energy in the exchange interac- 
tion. All wide binaries with orbital period P > 10 4 days are 
destroyed before our simulations begin. For each individual 
run, the mass of each BH is randomly selected with a distri- 
bution that reflects the results of BSR04, so that no two runs 
of any model contain the exact same BH population. 

The parameters used in all our simulations can be found in 
Table[T] For our simulations, aside from the exceptions noted 
in the table, we use self-consistent parameters determined by 
a King model, with Wq = 7, 9, and 1 1 . Given a total clus- 
ter mass, M c i, and core density, n c , we can calculate the one- 
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TABLE 1 
Simulation Parameters 





Structure 










n c 




ttb 


"I .core 


°"1.BH 


t'esc 


^halo.esc 


v halo 


Model Name 


(W ) 




Effective N 


Mm 


(pc - 


) 


(yr) 


(km s ) 


(km s ) 


(km s ) 


(km s ) 


(km s ) 


e5e5king7 .... 


7 


5xl0 5 


lx 


10 6 


512 


5x 


10 5 


1.5 x10 s 


14.1 


14.1 


55.6 


38.0 


40.7 


v2e5k7 


7 


5xl0 5 


lx 


to 6 


512 


5x 


10 5 


1.5xl0 8 


14.1 


7.0 


55.6 


38.0 


40.7 


v3e5k7 


7 


5xl0 5 


lx 


L0« 


512 


5x 


HP 


1.5 x10 s 


14.1 


4.7 


55.6 


38.0 


40.7 


v3e5k7ej54 a . . 


7 


5xl0 5 


lx 


to 6 


512 


5x 


HP 


1.5xl0 8 


14.1 


4.7 


55.6 


38.0 


40.7 


v3e5k7ej75 a . . 


7 


5xl0 5 


lx 


to 6 


512 


5x 


to 5 


1.5xl0 8 


14.1 


4.7 


55.6 


38.0 


40.7 


v5e5k7 


7 


5xl0 5 


lx 


10« 


512 


5x 


HP 


1.5xl0 8 


14.1 


2.8 


55.6 
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NOTE. — Starting from the top. the table is sorted according to the following parameters in their respective order: Wo, n L -, cti.bh, Mel, iVrjH- 

a These models include GW recoil, as described in i|3.4l Models v3e5k7ej54 and v3e5k7ej75 have maximum recoil velocities Vo = 54 and 75 km s _1 respectively. 
Models v2e55k9e6, v2e55k9e65, v2e55k9e7, and v2e55k9e8 correspond to V = 60, 65, 70, and 80 km s~' 
b These models include one primordial seed BH of mass IOOMq and 200Mq, respectively, as detailed in 13.31 
c The GMH model parameters correspond to the initial conditions used in GMH04. 



dimensional velocity dispersion, fTr core , the escape velocity 
from the center of the potential to the half-mass radius, v na io, 
the escape velocity of a BH in the core from the entire cluster, 
v eS c, and the escape velocity of a BH at the half-mass radius 
from the entire cluster, v na io,esc- We analyzed relatively mas- 
sive and dense clusters, precisely the types of clusters where 
BH growth would be expected. Specifically, we systemati- 
cally varied M z \ between 5 x 10 5 and 2 x 10 6 M Q , and n c be- 
tween 10 5 and 10 7 pc" 3 . Most of our simulations had 512 BHs 
(Nbh = 512), but in two simulations we looked at clusters with 
smaller and larger numbers of BHs. All of these parameters 
are listed in Tabled The escape velocities are used to de- 
termine whether the product of an interaction is to remain in 



the cluster as prescribed in 32.11 The velocity dispersion of 
the BHs, <tbh, can be related to the one-dimensional velocity 
dispersion simply as ctbh = 

The properties of the decoupled sub-cluster in relation to 
the properties of the cluster core are not so apparent. In our 
simulations we need to know the initial velocity dispersion of 
the BHs, <7bh- Numerical simulations suggest that the mean 
kinetic energy of the dynamically decoupled massive objects 
is only a few times larger than for other objects in the core 
(Gurk anet alJ20 04). When decoupling, the BHs will have an 
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FIG. 1. — Sub-cluster evolution. The data for each panel are binned and then averaged over several simulations in order to reduce the level of noise in the 
graph. The clusters reach equipartition when the total number of BHs in the core (solid) is < 40. Notice how few BHs are in the cluster halo (long-dashed) 
throughout the evolution. Only in v2e55kll, panel (c), do the numbers of BHs in the core and halo become comparable, but only when the cluster has already 
approximately reached equipartition. The more massive cluster v22e6e5k9, panel (a), does not reach equipartition in a Hubble time, suggesting massive clusters 
with low Wo could still contain significant numbers of single BHs in their cores. Each model has the same core density, » t =5x 10 5 pc~ 3 , and K = 4, except 
model v3e55k9, panel (d), which has K = 1.77. Model v22e6e5k9 is a massive W = 9 King model with M A = 2 X 10 s M . Model v2e55kl 1 is a W = 11 King 
model with M c j = 5 X 10 5 Mq. Models v2e55k9, panel (b), and v3e55k9 are the same as v2e55kll, except they have Wq = 9. 



effective velocity dispersion which we write as 

/ (m) y /2 

cr BH=[K— CT core , (11) 

where cr core is the velocity dispersion of the core, (Mbh) is 
the average BH mass, and K is the ratio of the mean kinetic 
energy of the BHs to that of the rest of the core, with K = 1 
corresponding to complete energy equipartition. Because the 
BHs are dynamically decoupled from the rest of the cluster, 
they will undergo their own independent evolution as a small 
cluster of stars, increasing their density and velocity disper- 
sion in the process. Therefore we look at sub-clusters with 
AT =.64, 1.77,4, and 16. 

3. CLUSTER DYNAMICS AND BH GROWTH 

In this section, we discuss our results for the evolution of 
the BH sub-cluster and we analyze the conditions under which 
successive mergers of BHs lead to the growth and retention of 
candidate IMBHs. 

In most of our simulations the BH sub-clusters reach 
equipartition (our standard for the complete disruption of the 
sub-cluster) in a few Gyr. In the most extreme cases the sub- 
cluster may dissipate in less than lOOMyr, as in the very dense 



model v3e5e7kl 1; or not at all, as in all models with K = 16. 
We also find there is a significant probability a BH with mass 
> 100M Q can form. For clusters that reach equipartition, 
massive BHs are most likely to form in clusters with high core 
densities and temperatures. This also leads to more growth 
for a given King profile. Also, King models with lower Wo are 
more likely to form candidate IMBHs. The growth of massive 
BHs is highly dependent on the maximum GW recoil veloc- 
ity, nearly stopping it completely even when it is only a few 
percent larger than the core escape speed. 

3.1. Fate of the BHs 

The fate of the BHs in the cluster is determined mainly by 
their characteristic interaction rate in the core. However, two 
main questions remain regarding how the space of initial clus- 
ter parameters is divided among the different possible out- 
comes: does the BH sub-cluster reach equipartition with the 
cluster stars within a Hubble time? Do the BHs experience 
enough successive mergers to form an IMBH? Because of the 
simplicity of our model and speed of our code, the trends as- 
sociated with varying one parameter in the simulation are ev- 
ident. 

Figure [2 shows the number of BHs located in the cluster 
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FIG. 2. — Mass distribution of largest remaining BH. The largest remaining BH is simply the most massive BH remaining in the cluster at equipartition. 
Panel (a) is the mass distribution from a few different cluster types, with the model name in the upper right-hand corner. Panel (b) is the mass distribution from 
cluster model v2e55k9, with varying GW recoil kick velocities, Vb, corresponding to, from the top of the figure down, models v2e55k9, v2e55k9e6, v2e55k9e65, 
and v2e55k9e7. 



core and halo as a function of time for a variety of different 
cluster types. We see that most BHs remain in the cluster core 
rather than in the halo of the cluster. A fraction < 3% of all 
the BH-BH mergers from cluster occur within the halo. The 
sudden drop in the binary fraction, as seen again in Figure^] 
can be attributed to strong binary-binary interactions. Even 
though this causes a low binary fraction for the majority of 
the dynamical evolution of the sub-cluster, approximately 10— 
15% of the BHs are ejected in binaries when the BH sub- 
cluster reaches equipartition in a Hubble time. 

For sub-clusters with low densities and high velocity dis- 
persions, the three-body binary formation rate can become so 
small that a newly formed binary is disrupted before the next 
binary formation. Despite this, our simulations suggest that 
there always exists at least one binary in the core that may be 
able to keep the sub-cluster from undergoing core collapse. 
This binary, among others, is created through three-body bi- 
nary formation, which is the underlying mechanism for the 
entire evolution of the sub-cluster. Since BHs are only ejected 
from the cluster after three- or four-body interactions, BH 
sub-clusters that do not produce binaries at a high rate do not 
dissolve within a Hubble time (see, e.g., models e5e5king7, 
e5king9, and v22e6e5k9). If the sub-clusters' parameters 
have not changed significantly over the evolution of the entire 
cluster, then a significant number of single BHs could exist in 
clusters similar to those considered in our simulations. 

3.2. Formation of an IMBH 

Our simulations uniquely allow us to follow the interactions 
of a realistic mix of single- and binary-BHs and monitor the 
growth of BHs through successive mergers. We find that in 
a given King model, clusters with greater core densities and 
larger masses have a greater probability of forming an IMBH. 
They usually grow to even larger masses. For a fixed mass 
cluster that fits a given King model profile larger core densi- 
ties result in a smaller half-mass relaxation time, t&. If this 
timescale is small enough, the cluster may actually undergo 
runaway growth through stellar collisions. This suggests that 



one is more likely to find an IMBH, whether formed through 
successive BH mergers or stellar collisions, in clusters with 
dense cores. 

However, clusters with lower degrees of central concentra- 
tion (smaller Wo, and higher <7 core ) seem to have more BH 
growth. Figure |2 shows the mass distributions of the largest 
BHs formed in a few different cluster types. We see that when 
clusters completely disrupt in less than a Hubble time, larger 
core temperatures directly correlate with more growth of a 
single massive BH. 

This is evident in the directly comparing models v2e55k9 
and v3e55k9 (where K = 4, and 1 .77 respectively). Although 
the difference in the average number of mergers in the cluster 
is small — v2e55k9 had about 14 mergers whereas v3e55k9 
had about 13 — the number and sizes of the large BHs formed 
are dramatically different. Model v2e55k9 had about twice 
the number of BHs as v3e55k9 with mass > 100M Q remain 
in the core at equipartition, and an average mass of the largest 
BH remaining in the core about 60% larger as well. In model 
v3e55k9, the velocity dispersion of the BHs, obh, is lower 
than v2e55k9. As is evident in equation (|9j, the interaction 
rate increases rapidly with lower velocity dispersions, and 
hence leads to a quicker dissipation of the BH sub-cluster. 
This is consistent with our understanding that three-body bi- 
nary formation drives the evolution to the eventual equiparti- 
tion of the system. 

In clusters where there are successive BH mergers, the first 
formation of BHs with mass > 100M Q occurs at ~ lOMyr 
after the sub-cluster formed, roughly independent of the clus- 
ter model. The probability that the cluster has a BH with mass 
> 1OOM is then proportional to the logarithm of time (until 
the cluster reaches equipartition). On such a short timescale, 
one may worry that the BHs could also be colliding with mas- 
sive stars, counter to our assumption of a "pure BH" system. 
However, as we discussed in § ^ f° r clusters in which run- 
away collisions of massive stars were avoided, the massive 
stars are not expected to enter the BH core before exploding 
as supernovae. 
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3.3. Proposed Mechanisms for BH Retention 

One explanation for why BHs would stay in the cluster and 
not be ejected by hardening and eventual recoil is the Kozai 
mechanism in st able hierarchical triples ( Miller & Ha milton 
2002a HWenll2 003). If the relative inclination between the or- 
bital planes in a given triple is large enough, the inner binary's 
eccentricity can be driven to a high value (formally to 1). The 
high eccentricity necessarily means a decrease in the merging 
time of the inner binary so that the inner binary will merge be- 
fore the triple's next likely interaction. We find that, typically, 
this has no significant effect on the formation of large BHs. 
Overall, mergers caused by the Kozai mechanism account for 
less than w 10% of the total mergers, and a negligible amount 
in most models. For example, the model with the largest per- 
centage of triple mergers, v3e5kl 1, exhibits almost no growth 
at all, with none of the 14 runs having a BH with mass greater 
than even 8OM . 

Another possibility that favors the retention of large 
BHs in clusters is the in troduction of a massive seed BH 
iMiller & Hamiltonl2002bl) . The mass spectrum of BHs given 
in BSR04 generally includes a BH of ~ 50M Q in each clus- 
ter. The presence of a BH of this mass does not always mean 
that this BH will remain in the cluster. In fact, the largest ini- 
tial B H is almost a l ways e jected from the cluster. Although 
iMiller & Hamilton! i2002bl) suggest that the introduction of a 
seed BH with mass 50M Q may be sufficient to cause signif- 
icant growth, we find that this is still not massive enough, in 
agreement with GMH04. 

In model v2e55k9-100 we introduce an initial seed BH of 
mass 100M Q ; in model v2e55k9-200 a 200M Q seed BH. 
We find that even BHs with these large masses can easily be 
ejected from the cluster. For example, in model v2e55k9- 
100, the seed BH was retained only 18% of the time. Even in 
model v2e55k9-200, with a seed BH of 200 M Q , the seed BH 
is again retained only 35% of the time. These probabilities are 
still lower than those in GMH04, where they found the BHs 
to be retained in a similar cluster ~ 40% and ~ 90% of the 
time, respectively. This discrepancy can likely be attributed 
to the mass distribution of BHs in our simulations. The av- 
erage mass of the BHs in our simulation is 50% higher than 
in GMH04, and thus there is an increased probability that the 
large seed BHs will be ejected. 

3.4. GW Recoil 

One key question regarding coalescing binary BHs is the 
magnitude of linear recoil caused by the asymmetry of the 
GW emitted by the binary. Our code allows us to prescribe a 
systemic recoil velocity for every BH-BH merger. Through 
this we can follow the effects of GW recoil in a cluster en- 
vironment. We are able to understand the possible effects of 
gravitational recoil by varying the maximum magnitude of the 
recoil velocity in the base model v2e55k9 (Favata et al .120041 
iBlanchet et alJ2005l) . 

As expected, the number of successive mergers has a strong 
dependence on the maximum recoil velocity; however, it 
seems that it has only a small effect on the overall dynam- 
ics of the rest of the BH sub-cluster. Increasing the recoil has 
almost no noticeable effect on the total number of mergers 
in the cluster, but can have significant consequences on the 
rate of BH-BH inspirals (see Even a maximum recoil 
velocity slightly larger than the escape velocity of the core 
(Yo = 1 .042 x v esc = 60 km s" 1 ) dramatically changes the num- 
ber of large BHs formed in this model, as seen in Figure [2] 




lo- 4 io-» 10-' 
Eccentricity 

FIG. 3. — Eccentricity distribution of merging BH binaries. For this 
log— log plot, we show the eccentricity distribution of all BH binary mergers 
for model v2e55k9. The distribution of eccentricities is almost entirely inde- 
pendent of the model used. The two frequencies of the GWs were chosen to 
show the expected eccentricity distribution of a binary as it enters the observ- 
able bands of both ground-based, 10 Hz, and space-based interferometers, 
~ 10~ 3 Hz. The low eccentricity of most binaries entering the ground-based 
interferometers detection band suggests almost no loss of detectable BH-BH 
binary signals if only circular templates are used for analysis. 



With this recoil velocity, the probability of forming a 100 M Q 
is cut in half. When we look at higher recoil velocities, the 
possibility of growing a large BH gets only smaller. 

This is not such a surprising result when one considers our 
simple prescription for modeling GW recoil. In mergers with 
a mass ratio, q = m.ilm.\ w .38, or when f(q) w f max , the co- 
alescing binary will generally be ejected when Vq is close to 
the core escape speed (see eq. J6)). For BHs which go through 
successive mergers, it is likely that before reaching masses 
> 100M Q the BH will have already been ejected. However, 
since this only affects BHs after they merge, the number of 
mergers remains relatively unaffected. 

Our simple model of GW recoil neglects some aspects of 
the process, which could alter results. In particular, we do not 
follow the evolution of the BHs' spins. Because of this, we 
must neglect the effects of spin on recoil, and therefore do not 
look at alternative paths to large BH retention. One possibil- 
ity is that if the BHs with the appropriate spins merge, some 
clusters could still retain larger BHs, even if the recoil veloc- 
ity of the BHs with no spin is much larger than the escape 
speed. On the other hand, spin breaks the symmetry of the 
binary and can lead to large reco il velocities even for equal- 
mass binaries (Favata et al.l2004l) . 

3.5. Properties of Merging BHs 

Inspiraling BHs with high eccentricity have very com- 
plex gravitational waveforms, and their additional free pa- 
rameters make computational searches even more expensive. 
Therefore we want to know the eccentricity of the merging 
BHs' orbits as they enter possible detection bands of cur- 
rent and future gravitational wave detectors. Many ground- 
based laser interferometers are currently in operation. These 
detectors, such as LIGO, are sensitive to GWs above » 
10 Hz, below which seismic noise dominates the noise curve 
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JCutler & FlanaganlfT99l . Also, in the planning stages is 
a space-based laser interferometer, LISA, which, with its 
longer baseline, is expect ed to be sensitive to GWs with fre- 
quenc ies ~ 10" 3 - 1 Hz dBender et alJ Il998t iBender & Hilsl 

Dynamical interactions, such as those in a cluster core, cou- 
pled with the strong dependence of merging time with eccen- 
tricity, suggest that many binaries will have highly eccentric 
orbits after their last strong interaction. Of course, GW emis- 
sion reduces the eccentricity of the orbit by circularizing the 
binaries (see eq. 0), and therefore it is often assumed that 
most binaries can be fitted with GW templates for zero ec- 
centricity. Figure|5]shows that there will be almost no loss in 
possible LIGO BH-BH binary sources with this assumption 
because, at such a high frequency, almost all binaries are cir- 
cular. However, the eccentricity distribution will likely matter 
for space-based detectors since the inspiraling binaries have 
not been entirely circularized by the GW emission. This is 
consistent with the previous study by GMH04. 

Another important factor in detecting inspirals is the chirp 
mass of a binary, 



(mim2) 3 / 5 
(mi +OT2) 1 / 5 



(12) 



which solely determines the overall magnitude of the GWs 
emitted by a coalescing circular binary. Because our sim- 
ulations allow for successive mergers of BHs, as shown in 
Figure |3 some mergers have chirp masses well above those 
expected if dynamics were not included. Because of our re- 
alistic initial distribution of BH masses, the chirp masses of 
most mergers is above the expected value for two 1OM BHs 
merging, 8.7M Q . 

Although, as we discussed in § 13.21 successive mergers of 
BH-BH binaries can lead to more inspirals of massive BHs, 
in most cases the chirp mass distribution at the end of the evo- 
lution of the cluster is not significantly different from early in 
the evolution. For one, about half of the mergers occur out- 
side of the cluster, which were ejected through dynamic inter- 
actions early in the cluster's history before significant growth 
had occurred. Also, only in a few cluster models is the prob- 
ability of BH growth near unity, in which the massive BHs 
dominate the mergers in the cluster. 

In Figure [5] we plot the average merger rate of model 
v2e55k9 as a function of time. The merger rate can clearly 
be broken into two regimes. The first occurs early in the clus- 
ter evolution, when the number of binaries has not been com- 
pletely depleted. The second regime occurs after the clus- 
ter is nearly all single BHs with only a few binaries. During 
this later time the merger rate of the cluster falls off inversely 
proportional to the age of the cluster, independent of cluster 
model. 

3.6. Comparison with Previous Studies 

Overall our results agree with those of previous 
direct jV-body simulations . In the simulations of 
iPortegies Zwart & McMillad feOOOl) . where N = 2048, 
4096 and Nbh = 20, 40, approximately 60% of the BHs where 
ejected as single BHs and 30% ejected as binary BHs. The 
lower values of binary ejection, which we found, can be 
explained by in-core mergers and the lower order interactions 
followed compared to A^-body simulations (the calculations 
of the authors were Newtonian only, and so did not allow 
for mergers via GW emission). As shown in Table [2] in 



most runs 70-85% of the BHs were ejected as singles and 
10-15% were ejected in binaries. 

The biggest divergence between our simulations 
and those conducted by Portegi es Zwart & McM illan 
(2000) is the energy distribu tion o f ejected BH binaries. 
Porteg ies Zwart & McMillan (2000) found that the ejected 
BH binaries had a binding energy distribution more or less 
uniform in the logarithm. In Figure [6] we show a typical 
binding energy distribution of the ejected BHs. Our simula- 
tions produce a distribution much more log-normal. Every 
model we analyze has a similar distribution of ejected BH 
binary binding energy, with values bet ween sa 10 3 - 10 5 kT , 
consistent with analytic considerations llKulkarni et alJl993h . 
This divergence can possibly be explained by the l ow number 
of BHs used in lPortegies Zwart & M cMillad ( 120001) . 

To see what effect the mass spectrum has on our simula- 
tions, and also to compare our data to the results of GMH04, 
we use three different models. Each simulation uses the same 
velocity dispersion and escape velocities as in GMH04 (see 
model GMH in Tabled, but starts with a different mass func- 
tion. GMHA and GMHC both have 10M Q and 15M equal- 
mass BHs respectively. GMHB has BHs with a mass distribu- 
tion consistent with our other simulations following BSR04. 
This distribu tion has a co rresponding average mass of about 
15M Q (see 32.1I & ^2.21 . As can be seen in Table|2j using 
equal-mass BHs in our simulations increases the number of 
binaries ejected by almost a factor of two. The models with 
equal-mass BHs, GMHA and GMHC, have about 20% of 
their BHs ejected as binaries, whereas GMHB has only about 
10%. FigureQshows how the mass spectrum we use causes 
the cluster to reach equipartition at an earlier time and changes 
the timescales of when mergers occur, compared with a clus- 
ter of equal-mass BHs. 

4. BINARY BHS AS SOURCES OF GRAVITATIONAL WAVES 

BH-BH binaries formed through dynamical interactions 
may be some of the best sources of GWs detectable by 
ground-based laser interferometers. Previous studies of de- 
tection rates have led to a large range of possible val- 
ues, with some of the greatest uncertainty coming from 
the dynamics of interacting BH binaries in massive clusters 
(Tutu kov & Yungelsor] Il993t IPortegies Zwart & McMillad 
2000). In this section we determine possible maximum de- 
tection rates of BH-BH binary inspirals from some globular 
cluster models. 

4.1. Calculation of the LIGO Detection Rate 

To calculate accurately the net detection rate, one should 
convolve self-consistent densities and birth rates of observed 
star clusters throughout the universe with the mergers rates in 
our simulations. Because there is great uncertainty in these 
distributions, we look at each of our cluster models, and de- 
termine how cluster densities and masses may affect the dis- 
tribution of BH-BH inspirals, leaving a detailed analysis for 
further study. 

Although current ground-based interferometers are not 
sensitive enough to detect mergers of inspiraling BH-BH 
mergers to any cosmologically significant distance, as these 
detectors become more sensitive they will be able to make 
detections to ever farther distances. It then becomes impor- 
tant that a consistent cosmological model is used to have a 
better understanding of detectable inspirals. In all of our 
calculation s we use the best-fi tting cosmological parameters 
found by Sperg eTet alJ (|2003) from combining WMAP data 
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FIG. 4. — Chirp mass of mergers versus time. This is a comparison of the two models v2e5kl 1 and v22e6e5k9, in panels (a) and (b) respectively. Plotted is the 
chirp mass versus time of all mergers of 46 random runs of model v2e5kl 1 and all 46 runs of v22e6e5k9. Model v2e5kl 1 is one of the least efficient clusters in 
producing large BHs and BH-BH binary mergers in general. Therefore, the distribution is most nearly that expected from the initial mass distribution of BSR04. 
Because of how quickly v2e5kll evolves (f cq ~ 200 Myr) almost all mergers in later times occur outside the cluster. In comparison, v22e6e5k9 is a massive 
cluster that does not reach equipartition before a Hubble time. There is still a significant fraction of BHs in the cluster at the end of the simulation, which allows 
for more growth, and also more massive BH mergers. 
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FIG. 5. — Merger rate vs. time. The solid curve is the average merger 
rate of model v2e55k9 as a function of time. The dotted line is a power- 
law oc time -1 . After ~ 1() 8 yr, the merger rate is inversely proportional to the 
age of the cluster. The evolution of the merger rates can be split into two 
phases. The first when the cluster is undergoing many binary interactions, 
and the second, when the binary fraction is depleted and nearly zero. These 
two phases of merger rates appear consistently in all cluster models. 
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FIG. 6. — Energy distribution of ejected BH binaries. Plotted is the proba- 
bility distribution of the energy of all BH-BH binaries ejected before equipar- 
tition in 117 runs of model v2e55k9. The energy is plotted in units of the 
mean kinetic energy kT, where 3/2kT is the mean stellar kinetic energy of 
the MS stars in the core of a cluster of this type. We find that all other models 
have a distribution very similar to the one shown above. 



with other measurements: Hq = 71 km s 1 Mpc Cl m = 0.27, 
ft 7 = 5 x 10" 5 , and tl A = .73. 

In our calculations, we assume that the globular cluster 
model was formed uniformly through the universe at a given 
cosmological time corresponding to redshift Zf orm . We then 
record each detectable merger into one of 100 bins each with 
time width St = fo/100, where to is the current age of the uni- 
verse, based on when the merger occurred. If d\ is the number 



of detections in bin i, we sum over the rate of each bin giving 
the final rate: 

100 , . 

* Zfo ™ = Y,j t Y po{D ' i ~ D ^ l){l+Zir ^ (13) 

where p is the current density of a given cluster model and z, 
is the redshift to bin i. With t e = fo - iSt, the proper distance to 
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FIG. 7. — A comparison of the models GMHA, GMHB, and GMHC. In panel (a), the number of BHs as a function of time is plotted. Each pair of lines is 
labeled above by its associated model that we use to compare our results with GMH04. Model GMHA, with only IOMq equal-mass BHs, on average reaches 
equipartition at f e q fa 5.2Gyr. Model GMHC, a cluster with only 15Mq equal-mass BHs, reaches equipartition at ( eq m 830Myr. Finally, GMHB, with a 
varied mass spectrum as in BSR04 and a corresponding average mass m 15Mq, has not only a significantly smaller number of ejected binaries, but also reaches 
equipartition at the earliest time, ( e q ~ 670 Myr. In panel (b), the number of mergers at a given time, denoted by redshift, is plotted. We assume that the clusters 
would currently be 13 Gyr old. The model with the varied mass, GMHB, not only dissociates more quickly, but has most of its mergers at an earlier time in the 
cluster's evolution, i.e., at a higher redshift. 



the edge of bin i is 



'<> dt 
0(0' 



(14) 



where a(t) is the scale factor of the Friedman-Robertson- 
Walker metric that satisfies the Einstein equation, and a(fo) = 
1. The factor (1 in equation Jl 3i comes from the cos- 

mological time dilation of the merger rate. 

The final detection rate is directly proportional to the den- 
sity of such clusters in the universe. For our calculations we 
assume that po = 1 Mpc" 3 , independent of cluster model, for 
ease of compari ng our results with other works. To put this 
in perspective, Porte gies Zwart & McMillan! (2000) found the 
number density of all globular clusters to be po sa 8 .4 h 3 Mpc -3 
based on rough fits to observations. Our analysis doesn't in- 
clude the full distribution of cluster parameters, but instead 
looks at each cluster individually. The Milky Way, for exam- 
ple, contains hundreds of globular clusters, with only a frac- 
tion o f clusters similar to those we look at in this study (lHarrisI 
1996). 

In order to be as precise as possible, we must determine 
which mergers could actually be detected by a given version 
of LIGO. To do this, we must look at the accumulated signal- 
to-noise ratio (SNR) of a given merger at redshift z m - Since, 
gravitational waveforms are invariant under redshift, we use 
the redshifted chirp mass of a given merger 



Mcimp = (l+Zm)M chiip , 
and its luminosity distance 

D h = (l+ Zm )D, 



(15) 



(16) 



to calculate the SNR of the merger, where D is the proper 
distance of the merger as calculated in equation dl4> . 

We look at all BH-BH binary mergers caused by interac- 
tions before the cluster reaches equipartition, and determine 



if it would be detected by a given ground-based GW interfer- 
ometer by comparing its SNR to that of an inspiraling neutron 
star-neutron star (NS-NS) binary at luminosity distance Z\,o- 
In practice, we determine the merger to be detectable if 



(S/A0(/off) = Dl.o 
(S/N)(D hfi ) D h 



M 



chirp 



M 



chirp. 



5/6 



/ *(/off) 
■S(/off,u) 



with 



f faB (/Q~ 7 / 3 , 

s(fos)= / c ; ^ df , 
Jo 



W) 



>1, (17) 



(18) 



where S n {f) is the detector's noise spectrum, A^ C Mrp,o is 
the effective chirp mass of the inspiraling NS-NS binary, 
and fnff.n is the cut-off frequency of the NS-NS merger 
( Cutler & Flanagan 1994). In this study, we assume that 
the mergers are isotropic and neglect the orientation of the 
merger relative to the detector. We use an analytic fit for 
the shape of Advanced LIGO's noise spectrum found in 
ICutler & FlanaganNl994 . 

„ ,,/. J oo f < 10Hz, ncn 

W )0C l(/o//') 4 + 2[l + (/V/o) 2 ] f > 10Hz, 

where fo = 70 Hz. We do not inclu de the coefficient as cal- 
culated by Cutler & Flan agan! i 19941) in this equation since it 
is completely canceled out in equation dl7t and rescaled by 
D L . Q . 

The final piece of equation d!7l > is the cut-off frequency of 
the merger, f s, after which the waveform of the GWs can 
not be accurately modeled. This frequency is generally re- 
garded to be the frequency of GWs at the binary's last circu- 
lar orbit (LCO), after which the BHs plunge into each other 
in a time less than the orbital p eriod. In their calculations, 
iKidder. Will. & Wisemanl lll993l) found that for two 1OM 
Schwarzschild BHs the frequency of the orbit at the LCO 
is 100 Hz, using what they called "hybrid" equations that 
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were valid through (post) 5//2 -Newtonian order for arbitrary 
masses. Because their calculations also showed that this is 
a lower limit of the orbital frequency of the binary for an ar- 
bitrary mass ratio, we choose to use this limit for calculating 
/off. For circular orbits, the GW frequency is twice the orbital 
frequency, therefore in our calculations we use 

/20M„\ / 1 \ 
/off «200(^j(— jHz (20) 

as the cut-off frequency of detectable GWs, where M = m\ + 
ni2 and z m is the redshift of the merger. The location of the 
LCO and its corresponding orbital frequency are far from well 
established. Nevertheless, numerical simulations and other 
analytic approximations have shown the orbital frequencies of 
equal-mass BHs to be only larger than the value we use here 
(Blanchet 2002J iGrandclement. Gourgoulhon. & Bonazzolal 
120021) . 

4.2. Results 

We see, in Figure|8] the expected detection rates of BH bi- 
nary mergers for cluster v2e55k9. Here we assume that all 
clusters of this model formed when Zf orm = 7. 84, o r 13Gyr 
ago, using the assumptions and equations of § 14.11 As can 
be seen in the figure, for the luminosity distances Advanced 
LIGO is expected to reach, the detection rate scales well by 
the power-law . This can be explained by the time evolu- 
tion of the merger rate of BHs. We find, for all simulations, 
that the merger rate scales as f" 1 after the disruption of the 
primordial binaries. For the distances of mergers Advanced 
LIGO is expected to be able to detect, the clusters are rela- 
tively old and hence, the rate changes very little. 

The detection rates and theoretical uncertainty of all mod- 
els can be found in the last two columns of Table [2] We find 
for model v2e55k9, for a version of LIGO capable of detect- 
ing NS-NS mergers at a luminosity distance, Dl.o ~ 190 Mpc, 
the net detection rate is w 2.7 yr" 1 . For v2e55k9e8, which 
has the same conditions as v2e55k9, but a recoil velocity 
Vb = 80 kms" 1 , the net detection rate is actually higher than 
v2e55k9 at as 4.1 yr -1 . One cause for this may seem to be 
the lower cut-off frequency of high mass binary mergers, but 
this is wrong. When analyzing the detection rate assuming a 
universal chirp mass for all mergers, the rate was still higher. 
It can actually be attributed to the higher merger rate later in 
the evolution of the cluster. This delayed merger rate may 
be a result of a few mechanisms. The dependence of merger 
rate in models v2e55k9-100 and v2e55k9-200 suggests that 
having more massive objects in the cluster results in a lower 
rate overall at the end of the evolution. Model v2e55k9-200, 
which has a 200 M Q seed BH, has a merger about half that 
of model v2e55k9. Another, slightly less significant mecha- 
nisms can still possibly be the masses of the merging binaries. 
The timescale for merger is longer when the masses of the ob- 
jects are less (see eqs. |4] and l5l). 

Within the level of theoretical uncertainty, the final detec- 
tion rate is roughly proportional to the initial number of BHs. 
This, of course, is not exact, and isn't expected to be. The 
timescale for cluster disruption is also proportional to the 
number of BHs, so it should be expected that there would be 
some variance in the merger rate at later times depending on 
the initial number of BHs. 

In order to get an estimated value of the actual Advanced 
LIGO detection rate, we must consider not only the number 
density of the massive globular clusters we look at here, but 
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FIG. 8. — Detection rate of BH-BH inspirals. The solid curve represents 
the expected detection rate of mergers from cluster v2e55k9 if it has a current 
density po = I Mpc~ 3 and formed at z\ = 7.84, or 13Gyr ago. The method 
and assumptions of our calculation are detailed in § |4] The dotted line is a 
power-law for R oc £>£ 3 . 



also the density of the low mass clusters. Our results suggest 
that the expected detection rate scales proportionally to the 
initial number of B Hs in the cluster. Considering th e num- 
ber of BHs that iPortegies Zwart & McMillan! J2000) quoted 
as being in the globular clusters they analyzed, we would ex- 
pect the actual rate to lie within the range of rates for the clus- 
ters reported in Table |2 rescaled to the number density of 
globular clusters in the universe (~ 3 Mpc" 3 ). 

5. SUMMARY AND DISCUSSION 

Our work is based on the reasonable assumption that, in a 
sufficiently large and dense star cluster, BHs created via stel- 
lar evolution concentrate in an inner core and effectively de- 
couple from the rest of the cluster following mass segrega- 
tion and the development of the Spitzer instability. Taking 
advantage of this decoupling, we have computed the dynami- 
cal evolution of the BHs in a highly simplified treatment of the 
stellar dynamics but covering a wide range of cluster models 
and repeating calculations for each model in order to obtain 
a complete statistical description of outcomes. Our assumed 
initial distributions of BH masses and binary parameters are 
based on the most recent population synthesis calculations for 
young stellar systems. In our simulations, we use a simple 
Monte Carlo method to follow the evolution of the BH subsys- 
tem in a fixed background cluster described by a King model. 
The BHs are allowed to interact only in the core, described by 
a constant density and velocity dispersion. All interactions in- 
volving binaries are computed exactly by direct (Fewbody) 
integrations but we implement three-body binary formation 
using a simple analytic rate formula (eq. (9))- Binary for- 
mation through dissipative two-body encounters is negligible 
in the systems we consider here (with velocity dispersions 
a < 100 kms" 1 ). We allow for Newtonian recoil of BHs into 
the halo, reintroducing them in the core following mass seg- 
regation. Between interactions, BH-BH binaries are evolved 
taking into account gravitational radiation and the possibility 
of a merger (eqs. |4] & 0). In some simulations, we attempt 
to model kicks due to gravitational radiation recoil in merging 
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binaries, parameterizing the large theoretical uncertainty in a 
simple analytic formula that depends on the binary mass ratio 
only (eq. Q). 

We present the results of our simulations in §|3]and0] and, 
in particular, we derive the probability of massive BH growth 
and retention in clusters (~ 20-80%). We also show that the 
Kozai mechanism in triples has almost no significant effect 
on the merger rate, BH growth, or retention, in contrast to 
previous suggestions in the literature. In addition, we derive a 
net rate of BH-BH binary mergers detectable by current and 
future ground-based GW interferometers. We find, under ex- 
tremely optimistic assumptions, that this can be up to a few 
tens of events per year, if we assume that the globular clusters 
likely formed over a time span Zf orm < 8. When in cluding re- 
coil from the gravitational "rocket" effect (Fitchett 1983) this 
rate is nearly doubled due to the increased merger rate from a 
population of less massive BHs. 

The simulations presented here improve upon previous 
studies in that we not only account for binary evolution due 
to GW emission between successive interactions of a fixed 
group of BHs, but we also use a realistic BH IMF based on 
the most recent population synthesis models (BSR04). Most 
analytic studies and numerical simulat ions so far assumed 
that all BHs were lOM^ fsee. e.g fMiller & Hamiltorl2002bl 
and Portegi es Zwart & McMillan 2000), or included just one 
large BH (GMH04). The mass spectrum from BSR04 gives 
a slightly higher average mass (ps 15M q ), and also includes 
at least one significantly large BH, which is almost always 
ejected early in our simulations. Even some of the largest 
BHs formed from mergers, with mass w 120M Q , are ejected 
from the cluster when they are formed early enough in the 
simulation to interact often with other large BHs. The distri- 
bution of the binary parameters also has a significant effect 
on the interactions. The low primordial binary fraction results 
in very few formed triples, and even fewer mergers in triples, 
whether enhanced by the Kozai mechanism or not. 

The results of our simulations indicate a greater likelihood 
of moderate growth in globular clusters than previous A^-body 
simulations have suggested possible, but they show varied re- 
sults compared with GMH04. Our wider BH mass function 
makes the ejection of BHs with masses as high as ~ 100M Q 
very likely. On the other hand it also encourages the formation 
of even larger BHs through successive mergers, and the chirp 
masses of merging binaries are larger than would be expected 
for a cluster with equal-mass ~ 10M© BHs. 

Despite the demonstrated growth, the probability for an 
IMBH of ~ 10 3 M Q to form directly via successive BH- 



BH mergers remains extremely small. However, we should 
not neglect the possibility of significant further growth of 
the final remaining BH through stellar collisions after com- 
plete evaporation of the BH sub-clust er. Results from 
Baumgardt. Makino. & Ebisuzaki ( 2004a) suggest that even a 
moderately massive ~ 200M Q BH could grow into a 10 3 M Q 
IMBH after only a few Gyr, well within the current ages of 
globular clusters. 

Initial conditions play a crucial role in determining the 
probability of IMBH growth. In general, massive, dense clus- 
ters with high core temperatures have the greatest likelihood 
of BH growth. Eccentricity growth through the Kozai mech- 
anism, although it makes it possible for BHs to merge with 
very little (Newtonian) recoil, does not occur often enough 
to affect IMBH formation. Note, however, that our basic as- 
sumptions break down at late times, when the growth of BHs 
is most significant. Eventually the number of remaining BHs 
becomes small enough that the constant core conditions as- 
sumed in our simulations are no longer justified. In particular, 
the BH subsystem will start interacting with lower-mass stars 
at a significant rate and it will then likely return to energy 
equipartition. This is why we emphasize the values computed 
by our code at equipart ition , as opposed to the values at the 
end of the evolution (§ 11.21 . Beyond the return to equipar- 
tition, the evolution of the remaining BHs will be strongly 
coupled to the overall evolution of the whole star cluster. 

Our simulations are far from complete, yet the results show 
great promise for future work. For one, we can expect that the 
most massive clusters will have the greatest likelihood of BH 
growth. To enter the next regime of even more massive clus- 
ters, such as galactic nuclei, one must account for bin ary for- 
mation from gravitational bremsstrahlung (lLedll993l) . Also, 
because of the simplifications we make, it is also possible to 
explore an even larger parameter space, especially of the more 
numerous and less massive clusters. By convolving these data 
with the cluster formation history of the universe, one could 
determine a BH-BH merger detection rate to an ever more 
accurate degree. 
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NOTE. — Results of all simulations. Starting from the left the columns are the model name, number of total runs made, the average mass of the largest remaining BH at equipartition, the standard deviation 
of the largest mass, the total number of successive mergers the average largest remaining BH went through before equipartition, the largest BH formed that remained in the cluster, the fraction of simulations 
with a BH off mass > IOOMq present at equipartition, the fraction of BHs ejected as singles (/e.sin)> the fraction of BHs ejected in binaries ( fc.bin), the log of time at which the cluster reached equipartition(r cq ), 
the fraction of mergers in the cluster that happened in triples (N lT \p/N c \ uslCT ), the fraction of total mergers that occurred in the cluster (A' C | us[i;r /A^ mer g C ), the average total number of mergers for a single simulation 
(Amcrge), the expected Advanced LIGO rate of detection if the cluster formed when ^f orm = 7.8 (^7.g), and the expected Advanced LIGO detection rate if the cluster formed when zr orm = 1 (R\). For our calculation 
of the expected detectable merger rate, we assume that the given cluster model has a current uniform number density po = 1 Mpc~ 3 and that the interferometer is capable of detecting a NS-NS inspiral up to 
Dl,o = 190 Mpc (see ^4]for our detailed calculations and assumptions). 



a In these simulations the cluster never reaches equipartition before a Hubble time. All numbers are calculated at the end of the simulation at log 1( 



: 10.13. 



